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ABSTRACT 

With the detection of extrasolar moons (exomoons) on the horizon, it is important to consider 
their potential for habitabihty. If we consider the circumstellar Habitable Zone (HZ, often 
described in terms of planet semi-major axis and orbital eccentricity), we can ask, "How does 
the HZ for an Earth-like exomoon differ from the HZ for an Earth-like exoplanet?" 

For the first time, we use ID latitudinal energy balance modelling to address this question. 
The model places an Earth-like exomoon in orbit around a Jupiter mass planet, which in turn 
orbits a Sun-like star. The exomoon's surface is decomposed into latitudinal strips, and the 
temperature of each strip is evolved under the action of stellar insolation, atmospheric cooling, 
heat diffusion, eclipses of the star by the planet, and tidal heating. We use this model to carry 
out two separate investigations. 

In the first, four test cases are run to investigate in detail the dependence of the exomoon 
climate on orbital direction, orbital inclination, and on the frequency of stellar eclipse by the 
host planet. We find that lunar orbits which are retrograde to the planetary orbit exhibit greater 
climate variations than prograde orbits, with global mean temperatures around O.IK warmer 
due to the geometry of eclipses. If eclipses become frequent relative to the atmospheric ther- 
mal inertia timescale, climate oscillations become extremely small. 

In the second investigation, we carry out an extensive parameter study, running the model 
many times to study the habitability of the exomoon in the 4-dimensional space composed of 
the planet semi-major axis and eccentricity, and the moon semi-major axis and eccentricity. 
We find that for zero moon eccentricity, frequent eclipses allow the moon to remain habitable 
in regions of high planet eccentricity, but tidal heating severely constrains habitability in the 
limit of high moon eccentricity, making the habitable zone a sensitive function of moon semi- 
major axis. 
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1 INTRODUCTION 

The increasing rate of discovery of extrasolar planets (exoplanets) 
has given astronomers licence to consider their potential as habi- 
tats for extraterrestrial life. To date, more than 860 exoplanets have 
been confirmeqj, with over 2,300 exoplan et candidates identified 
by the Kepler mission jBatalha et alj|2013h . Of the confirmed ex- 
oplanets, several are in the circumstellar Habitable Zone (HZ), an 
annular region surrounding the star where planets with orbits inside 
it may be expected to possess liquid water on their surface. 

The HZ is derived for planets of terrestrial mass and com- 
position - the iimer edge is determined by water loss (via hy- 
drogen escape and photolysis), and the outer edge is deter- 
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m ined b y runaway glaciation due to a build up of CO2 clouds 
(Kasting, Whitmire & Reynolds 1993). In the case of single stars, 
the HZ is a circular annulus, and hence terrestrial planets on circu- 
lar orbits inside the annulus are expected to be habitable. Even if 
the planet's orbit is not circular, it may still be habitable depend- 
ing on the average received flux over its orbit, which can vary sig- 
nificantly if the orbit i s highly elliptical (Williams & Pollard 2002|; 
I Kane & Gelincl2012a la). If the star is a multiple system, the geom- 
etry of the HZ annulus becomes time-dependent, and this can have 
important consequences for the c limates of planets in t he system 
( lForganl2012l : lEggl et al.ll2012bllat iKane & Hinkeill2013l^ . 

To date, all exoplanets that reside within their local single- 
star HZ have radii 1.5 to 2 times larger than that of the Earth, 
and as s uch may be mini- Neptunes than super-Earths (see for 
example I Sames et alj 120091) . Even if such planets were terres- 
trial in nature, they are likely to have significantly higher man- 
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tie viscosity, which suppresses convection and as a result has 
implications for heat transfer and the planet's ability to degas 
JStamenkovic. Spohn & Breuedl2010l. iLiu & Zhon j|2010l) . 

This being the case, the current crop of exoplanets in the HZ 
should not be ruled out as a source of habitable niches. Extrasolar 
moons (exomoons) belonging to these planets may be terrestrial in 
natur e, and as a result may be habitable Jwilliams. Kasting & Wadd 
1 19970 . This is true even in our own Solar System: outside the Earth, 
the best chances for life in the Solar System are believed to be on 
tidally heated icy moons which are thought to pos sess liquid sub- 
surface oceans, such a s Europa (IMelosh et al .1120041) and Enceladus 
dParkinson et al.l2007t) . Titan, the largest moon in the Solar System, 
has a thick atmosphere, with evi dence of precipita tion and lakes 
constituting a hydrological cycle JStofan et al.l2007D . 

Exomoons are yet to be detected, but Earth-mass ex- 
omoons are in principle detectabl e indirectly with Kepler 
JKipping. Fossev & Campanella 120091) . Studying an exoplanet's 
transit timing variations and transit duration variations (TTV and 
TDV respectively) allows observational constraints on the moon's 
mass, and its orbital semi-major axis around the host exoplanet. Ad- 
ditionally, auxiliary transits of the moon and planet-moon mutual 
events lead to unique transit features revea ling the moon's radius 
and other orbital parameters (lKippind201lh . Such efforts to realise 
a detection are under way as part of the H unt for Exomoons with 
Kepler (HEK) project JKipping et al]l2012l') . 

A large component of habitability studies is the study of the 
energy balance associated with the celestial body in question. The 
various sources and sinks of heat are catalogued and compared, and 
the system is evolved over time to understand their interplay, either 
by analytical calculation or numerical simulation. In the case of ex- 
omoons, the list of sources is longer than that of exoplanets. Along- 
side direct stellar insolation, exomoons receive starlight reflected 
from the planet (as well as the planet's own radiation), and they are 
tidally heated. Also, the moon is more prone to stellar eclipses as a 
result of the host planet, providing an extra (effective) energy sink. 

The study of habitable exomoons has grown apace in recen t 
years: early calculations by [Reynolds. McKav & Kasting (19871) 
demonstrated that Europa is a viable niche for marine life such as 
that found under Antartic ice, and proposed a tidally heated circum- 
planetary habitable zone. In the ea rly phase of exoplanet detection, 
Iwilliams. Kasting & Wadd 119971 ) suggested that two of the then 9 
known giant exoplanets could host habitable moons. They raised 
the three key obstacles to exomoon habitability, namely: 

(i) Tidal locking of the moo n with respec t to the planet, resulting 
in extreme weather conditions ( iJoshill 19970 . 

(ii) Increased EUV and X-Ray radiation from the host 
planet's magnetosp here leading to loss of the moon's atmosphere 



planet s magnetosp nere 1 
( lKalteneggedl20 0Cf), and 



(iii) A poorly proportioned distribution of volatile abundances, 
due to the d iffering environments of planet and moon formation. 
iKalteneggen yOOOf) notes that Solar System moons may be over- 
abundant with water to be habitable, if placed in the inner Solar 
System, as they would lack the large areas of land required to sus- 
tain a carbonate-silicate cycle. Without this, the body would lack a 
crucial regulatory system to ameliorate the effect of the star's grow- 
ing luminosity as it evolves along the Main Sequence. 

IScharn ( 120061) investigated the then-known exoplanet popula- 
tion for the potential to host habitable exomoons. Estimating that 
almost 30% of the total population could host small moons with 
icy mantles, he goes on to estimate that tidal heating could increase 
the habitable zone's extent in planet semi-major axis by a factor of 



around 2 (although he notes that the tidal heating required can be- 
come several orders of magnitude larger than that experienced by 
lo, o ne of the more ex t reme c ases in our Solar System). 

iHeller & BamesI 1 120131) points out that exomoons may in 
some cases provide better niches than exoplanets, for several rea- 
sons: 

(i) Exoplanet s in the H Z of low mass stars (e.g. 
[Dressin g & Ch arbormeaul 1201 3b are likely to be tidally locked, 
whereas moons of such exoplanet s are likely to be tidally 
lock ed to the planet, not the star dKaltenegger & SelsisI l20ld : 
Sasa ki. Barnes & O'Brienll2012r) . This reduces the strong fluctua- 
tion of surface insolation (and potential atmospheric collapse) that 
a moon would suffer if it was tidally locked to the star (Joshi 19971 : 
iKite, Gaidos & Manga 1201 ih . This may also provide energy to 
drive the moon's internal dynamo, and generate a magnetosphere 
to guard against atmospheric loss. 

(ii) Exoplanet obliquity will erode due to tidal evolution, as 
demonstrated by the axial orientations of Venus and Mercury. This 
erosion is resisted more easily by more massive planets, if its 
orbital semi major axis is not too low and the planet does not 
posse s s relatively massive neighbours dHeller. Leconte & BamesI 
1 20 III : iSasaki. Barnes & O'Brien! 120121) . As a result, exomoons 
tidally locked to their host planet can retain the obliquity of their 
host planet and not suffer such strong erosion. 

In their work, JHeller & BamesI ( |20I3|) calculate the time de- 
pendent flux upon a moon due to the star's insolation, reflected 
starlight and thermal emission from the planet's surface, and add 
to this the tidal heating received. They use this total effective flux 
to calculate a parameter space of habitability in the moon's orbit 
about the host planet, and the host planet's mass, which they in- 
vestigate for Kepler-22b and KOI211.01, both of which reside in 
the circumstellar HZ. By averaging over the moon's orbital period, 
they produce a corrective term for the circumstellar HZ depend- 
ing on the planet's albedo and the moon's orbital semi-major axis 
relative to the planet, as well as surface temperature maps for the 
moons. However, while these calculations are among the most so- 
phisticated to date, they do not allow for the transfer of heat due to 
an atmosphere. This is clearly a crucial component of any climate 
model, particularly in the case of exomoon systems, which present 
many opportunities to set up strong temperature gradients that will 
generate heat transport. 

In this paper, we implement a ID latitudinal energy balance 
model (LEBM) to study the climate of an Earth-like exomoon in 
orbit around a Jupiter mass planet, which in turn orbits a Sun-like 
star. Our use of a LEBM allows us to go beyond analytical calcu- 
lations of detailed energy balance by allowing the moon to transfer 
heat across its surface (via a simple diffusion approximation) as 
well as adding some of the sources and sinks of energy described 
above. 

Rather than investigate the circumstellar HZ separately from 
the circumplanetary habitable region, we choose instead to investi- 
gate them together as a combined four dimensional manifold, con- 
stituted by the planet's semi-major axis and eccentricity relative to 
the star (a^, e^) and the moon's semi-major axis and eccentricity 
relative to the planet {am, e,„). We wish to investigate how the extra 
dimensions of this manifold produces features unique to exomoon 
habitability, and how other factors, such as the moon's orbital di- 
rection, affect the resulting climate. 

We investigate this problem in two parts: in the first, we sim- 
ulate four simple test cases to study in detail the effect of orbital 
dynamics on the resulting climate. In the second, we carry out a 
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large parameter study of simulations to investigate the four dimen- 
sional habitable manifold. 

In section [2] we describe the LEBM and its implementation; 
in section[3]we display results from our investigations; in section|4] 
we discuss the results, and in section|5]we summarise the work. 



2 LATITUDINAL ENERGY BALANCE MODELLING 

2.1 Simulation Setup 

In all simulations, we fix the star mass at M* = lAf©, the 
mass of the host planet at Mp = lA/jup, and the mass of the 
moon at IA/q. This system has been demonstrated to be dynam- 
ically stable on timescal es comparable to the Solar System lifetime 
JBames & OBrienl2002l) . 

The planet's orbit around the star is specified by its semi-major 
axis Op and its eccentricity Cp, and the moon's orbit around the 
planet is specified by its semi-major axis a™, eccentricity Cm and 
inclination relative to the planet's orbital plane im (see Figure [TJ. 
The orbital longitudes of the planet and moon are defined such that 
(pp = 4''m = corresponds to the x-axis. Note that as the moon 
mass is much less than the planet mass, we assume am is equal to 
the distance of the moon from the planet, as this is approximately 
equal to the distance from the moon to the barycentre of the planet- 
moon system. 

Comparing this to the more traditional setup for exoplanet 
LEBMs, there are two principal differences: 

(i) The epicyclic motion of the moon relative to the star (as de- 
scribed by the distance between the star and the moon, r«m), and 
(ii) The potential for stellar eclipses by the planet. 

A third difference is the potential for an extra energy source in the 
fo rm of tidal hea ting, which we approximate using the expressions 
in lScharfl | |2006[). While ou r calculations are of a lower dimension 
than those of Heller & Barnes! (|2013|)By allowing for diffusion o f 
heat across latitudes, we also differ from lHeller & BamesI (120130 . 
We describe the equations that constitute the LEBM in the follow- 
ing section. 



2.2 A One Dimensional Latitudinal Energy Balance Model 
witii Tidal Heating 

In this paper, we adapt the one dimensional LEBM applied to plan- 
ets to function for moons with an atmosphere of similar composi- 
tion to the Earth. LEBMs solve the following diffusion equation: 



C 



dT 
'dt 



d_ 
dx 



D{1 



2,dT 



= S[1-A{T)]^I{T), 



(1) 



where T = T{x, t) is the temperature at time t,x — sin A, and A is 
the latitude (between —90° and 90°). This equation is evolved with 
the boundary condition ^ = at the poles. The (1 — x'^) term is 
a geometric factor, arising from solving the diffusion equation in 
spherical geometry. 

C is the eifective heat capacity of the atmosphere, _D is a diffu- 
sion coefficient that determines the efficiency of heat redistribution 
across latitudes, S is the insolation received from the star, I is the 
atmospheric infrared cooling and A is the moon's albedo. In the 
above equation, C, S, I and A are functions of x (either explicitly, 
as S is, or implicitly through T{x, t)). 

D is defined such that an Earth-like planet at 1 AU around 
a star of IMq, with diurnal period of I day will reproduce 



the average tem perature profile measured on Earth, see e.g. 
ISpiegel. Menou & Scharf ( 2008) . Planets that rotate rapidly expe- 
rien ce inhibited l atitudinal h eat transport, due to Coriolis eff ects 
(see lFarrell|[l990h . We follow ISpiegel, Menou & Scharfl (l2008b bv 
scaling D according to: 



D = 5.394 X 10^ 



^d 

^d,® 



(2) 



where ujd is the rotational angular velocity of the planet, and uid,s, 
is the rotational angular velocity of the Earth. This is a simplified 
expression, which does not describe the full effects of rotation, but 
allows for rapid computation without severely compromising the 
model's accuracy. 

The diffusion equation is solved using a simple explicit for- 
ward time, ce ntre space finite difference algorithm (as was done in 
lForgai]||20I2h . A global timestep was adopted, with constraint 



St < 



{Axf C 
2D(l-x2; 



(3) 



The parameters are diumally averaged, i.e. a key assumption 
of the model is that the moons rotate sufficiently quickly relative to 
their orbital period. This is clearly not applicable for certain orbital 
param eters, as tida l locking will play a significant role at low values 
of a„ jJoshilll997h . 

The atmospheric heat capacity depends on what fraction of 
the planet's surface is ocean, f ocean, what fraction is land fiand = 
1.0 — f ocean, and what fraction of the ocean is frozen fice'- 



^ — Jland^land r Jocean \\^ Jice J^ocean ~r Jice^ice) • 

The heat capacities of land, ocean and ice covered areas are 
Ciand = 5.25 X 10^ erg cm-^ Yr\ 



(9.2Czand 
0.0 



263 K < T < 273 K 
T < 263 K 

T > 273 K. 



The infrared cooling function is 

UT) = ^^^ 

^ ' 1 + 0.75tir{T)' 

where the optical depth of the atmosphere 



riR{T) = 0.79 



273 K 



The albedo function is 
A{T) = 0.525- 0.245 tanh 



T - 268 K 



5K 



(4) 

(5) 
(6) 

(V) 

(8) 

(9) 
(10) 



This produces a rapid shift from low albedo to high albedo as the 
temperature drops below the freezing point of water, producing 
highly reflective ice sheets. It is this transition that makes the outer 
habitable zone extremely sensitive to changes to various orbital 
and planetary parameters, and makes LEBMs an important tool in 
studying short-term temporal evolution of planetary climates. 

The insolation flux 5* is a function of both season and latitude. 
At any instant, the bolometric flux received at a given latitude at an 
orbital distance r is 



S — qo cos Z 



lAU 



(11) 
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Figure 1. The setup for the latitudinal energy balance model (LEBM), shown in the x-y and x-z planes respectively. We assume the moon's mass is negligible 
compared to the planet mass, and as such the moon orbits the planet, not the barycentre of the planet-moon system. 



where go is the bolometric flux received from the star at a distance 
of 1 AU, and Z is the zenith angle: 

4 

(12) 



go = 1.36 X 10'' ( ^ I ergs - cm 

cos Z = ^ = sin A sin 5 + cos A cos 5 cos h. 



(13) 



We have assumed a simple main sequence scaling for the luminos- 
ity. 5 is the solar declination, and h is the solar hour angle. The 
solar declination is calculated from the obliquity So as: 



sin 5 



sin(5o cos( 



(14) 



where (f)trn is the current orbital longitude of the moon relative 
to the star, (fjperi.m is the longitude of periastron, and (fja is the 
longitude of winter solstice, relative to the longitude of periastron. 



We set ( 



= <j)a = for simplicity. 



We must diurnally average the solar flux: 



<?o^- 



(15) 



This means we must first integrate fi over the sunlit part of the 
day, i.e. h — [~H,+H], where H is the radian half-day length 
at a given latitude. Multiplying by the factor H/tt {&& H = -k if 
a latitude is illuminated for a full rotation) gives the total diurnal 
insolation as 

S — qo [ — ] fi ~ — (H sin A sin 5 + cos A cos 5 sin H) . (16) 

\ TT / TT 

The radian half day length is calculated as 

cos iy = — tan A tan (5. (17) 

In the interest of computational expediency, we make a simple ap- 
proximation for tidal heating, by firstly assuming the tidal hea ting 
per unit area is jPeale, Cassen & Revnoldslll98(]0scharfl2006l) : 



er = 



21 pl,R 



5/2g2 



38 



TQ 



5/2 



(18) 



where T is the moon's elastic rigidity (which we assume to be uni- 
form throughout the body), Rm is the moon's radius, pm is the 
moon's density, Mp is the planet mass, am and e^ are the moon's 
orbital semi-major axis and eccentricity (relative to the planet), and 
Q is the moon's tidal dissipation parameter. We assume terrestrial 
values for these parameters, hence Q — 100, F — 10^^ dyne cm~^ 
(appropriate for silicate rock) and pm. = 5 g cm"''. 

To calculate how this heating is distributed latitudinally, we as- 
sume the heating rate is at its maximum at the subplanetary point. 



We reuse equation ( I16t . substituting go for ex, S for the appro- 
priate S of the moon relative to the planet (which in this case is 
equal to zero as we assume So, i.e the relative obliquity between 
the moon and the planet, is zero), and fixing H — tt (i.e. tidal heat- 
ing occurs throughout the diurnal period of the moon). This is very 
much an approximation - indeed, the multi-dimensional nature of 
tidal heating prohibits latitudinal models from improving much on 
approximations such as this - habitability studies typically assume 
tidal heating is uniform ly distributed throughout the body (see e.g. 
lHeller&BMnesll2013l) . 

We calculate habitabili t y ind ices in the same manner as 
ISpiegel. Menou & ScharSfcoOSi) . The habitability function ^ ifl 



?(A,t) 



1 273K<r(A,i) <373K 

otherwise. 



(19) 



We then average this over latitude to calculate the fraction of hab- 
itable surface at any timestep: 



C(i) = ^r C(A,i)cosAdA. 

^ J-tt/2 



(20) 



Each simulation is allowed to evolve until it reaches a steady or 
quasi-steady state, and the final ten years of climate data are used 
to produce a time-averaged value of S,{t), 5, and the sample stan- 
dard deviation, a^. We use these two parameters to classify each 
simulations as follows: 

(i) Habitable Moons - these moons possess a time-averaged ^ > 
0.1, and a^ < 0.1^, i.e. the fluctuation in habitable surface is less 
than 10% of the mean. 

(ii) Hot Moons - these moons have temperatures above 373 K 
across all seasons, and are therefore conventionally uninhabitable 

(e<o.i). 

(iii) Snowball Moons - these moons have undergone a snowball 
transition to a state where the entire moon is frozen, and are there- 
fore conventionally uninhabitable (^ < 0.1). 

(iv) Transient Moons - these moons possess a time-averaged 
(, > 0.1, and a( > 0.1^, i.e. the fluctuation in habitable surface 
is greater than 10% of the mean. 



^ This function is often labelled r; - we choose g to avoid confusion with 
the use of r] as the frequency of Earth-like/habitable planets used in the 
exoplanets/SETl communities 
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3 RESULTS 

Exomoon systems present a parameter space of high dimension. We 
therefore take two approaches to exploring subsets of this space. 
The first approach fixes the planet's orbit, and considers four dif- 
ferent sets of moon orbital parameters as test cases. 

The second is a larger survey of parameter space, systemati- 
cally varying ap,ep,am,em (with the other orbital parameters held 
fixed). This allows us to map out a four-dimensional manifold to 
represent the habitable zone, a subset of the true, higher dimen- 
sional habitable zone manifold (which will depend on extra param- 
eters such as the obliquity of the moon and its chemical composi- 
tion, which we do not vary). 



3.1 Four Test Cases 

We fix the orbit of the planet with Up = 0.9 AU with e^ — 0, 
and consider four simple test cases to probe the effects of epicylic 
motion and eclipses on exomoon habitability. 

(i) A prograde circular orbit (i.e. where the orbital angular mo- 
menta of the planet and moon are aligned), with a™ ~ 0.023AU = 
iSlRp (approximately 0.3 Hill Radii), with im = 90° 

(ii) As 1., but with a retrograde orbit (where the orbital angular 
momenta are anti-aligned, and im ~ 270°). 

(iii) As 1., but with a semi-major axis reduced by a factor of 5: 



„, = 0.0046AU = 96Rp 
(iv) As 1., but with im -- 



(approximately 0.06 Hill Radii). 
= 0° (i.e. a polar orbit). 



These parameter choices are not entirely physically motivated - for 
example, it is unlikely that the or bit described in case 3 w ould re- 
main stable on long time scales (Weidn er & HomduOlOl) , not to 
mention that synchronous rotation would become important to the 
moon's climate. We choose these cases to illustrate more clearly the 
dynamical effects of changing the orbital parameters of the moon. 
We select ap — 0.9 AU to guarantee that all four test cases are 
habitable (see following section). 

All four moons are entirely habitable (^ — 1). This is easily 
established by studying the mean, minimum and maximum surface 
temperatures (Figure |2). The minimum temperature does not de- 
crease beneath 290 K throughout the orbit of the moon and the 
planet. The differences between prograde and retrograde moons 
are immediately obvious (top left and top right panel respectively). 
While the mean temperatures appear to be similar, the maximum 
and minimum temperatures fluctuate much more strongly in the ret- 
rograde case. Case 3 is very similar to case 1, and case 4 also has 
a similar mean temperature (with a more pronounced oscillation in 
the maximum and minimum temperatures). 

Why should we see such a large change by changing some- 
thing as simple as the direction of orbit, and yet see little change by 
reducing the moon's orbital distance by a factor of 5? The answer 
lies in considering the mean temperature more carefully (Figure[3]l. 
The mean temperature fluctuates by < 0.1 K, in all cases, but it is 
the nature of the fluctuations that are telling. The retrograde orbit 
is consistently hotter than the prograde orbit, and the frequency of 
temperature fluctuations is also higher. Case 4, where the exomoon 
orbits perpendicular to the plane, exhibits a superposition of fre- 
quencies, and case 3, with low am, exhibits essentially no variation 
at all. 

We can find an explanation for these phenomena if we con- 
sider the orbital dynamics. Let us consider first the difference be- 
tween the prograde and retrograde cases. 




Figure 4. The precession of the moon's longitude of periastron. The sun and 
planet are marked in yellow and red respectively, with the moon's longitude 
of periastron indicated by the blue dot. 



The planet's orbit is circular, and as such has no apastron or 
periastron. The moon's orbit is also circular relative to the planet, 
and hence has no apojove or perijove. But, relative to the star it 
undergoes epicyclic motion, and hence the moon has well-defined 
periastron and apastron. Figure |4]shows the location of the moon's 
periastron over the course of the planet's orbit. 

The longitude of moon periastron, ^pcr, precesses with a pe- 
riod equal to the planet's orbital period, Ppi. We define the x-axis as 
corresponding to = 0, and therefore initially, (f>pcr = vr radians. 
So, 



0pcr(i) ~ TV + (fipit mod 27r. 



(21) 



If the planet were static ((^pi — 0) then the longitude of periastron 
would not precess, and the time taken for the moon to move from 
apastron to periastron would be 



To = 



(22) 



where (pm — 27r/P,„ = n,„, and Pm is the moon's orbital period. 
However, ^pi 7^ 0, and therefore the longitude of periastron moves 
during this interval. If the moon's orbit is prograde, the angular dis- 
tance between apastron and periastron increases, and the timescale 
becomes 



rp 



TT + (pplTp 



If the orbit is retrograde, the angular distance decreases, and 



TR 



71" — (pplTR 



Rearranging the above equations gives 

TT 



rp = 

and 

rp = 



»m — <?'pl 



(23) 



(24) 



(25) 



(26) 



respectively. Therefore, the ratio of one to the other is 
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Figure 2. The temperature evolution of the four test cases, over 2 years, approximately 18 moon orbital periods (560 for case 3), and approximately 2.35 
planetary orbital periods). The top left panel shows case 1, the top right panel shows case 2, the bottom left panel shows case 3, and the bottom right panel 
shows case 4. The solid lines indicate the mean temperature, the dashed Hnes indicate the global minimum temperature, and the dotted lines indicate the global 
maximum temperature. 
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Figure 3. The mean temperature of the four test cases, over 2 years, approximately 1 8 moon orbital periods (560 for case 3), and approximately 2.35 planetary 
orbital periods). The top left panel shows case 1, the top right panel shows case 2, the bottom left panel shows case 3, and the bottom right panel shows case 
4. 
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and as such the moon's periastron and apastron distances become a 
function of the orbital longitude of the planet. 

At i = 0, the periastron and apastron distances are equal, and 
the distance between the star and the moon r,m will not change dur- 
ing its orbit. After 0.25 planetary orbital periods have elapsed, the 
periastron and apastron distances are now equal to that in the other 
3 cases. Similarly, at i = 0.5 planetary orbital periods the perias- 
tron and apastron distances are the same, and at t = 0.75 planetary 
orbital periods they reach their maximum values once more. This 
provides the amplitude modulation seen in the bottom right panel 
of Figure[3] Note the frequency of the oscillation is the same as that 
of case 1, as they are both prograde orbits. 



Figure 5. Arrangement of Lp and L^ , the planet and moon orbital angu- 
lar momentum vectors respectively, in test case 3, a polar orbit. The figure 
shows the arrangement at time t = 



TP 



4>Yn — 4>pl 
4>in + 4>p\ 



(27) 



The orbital peri ods of both bodies are related by (Appendix B of 
lKi ppin^l20()9J) : 



pi 



(28) 



where x is o-m normalised by the moon's Hill Radius. This gives 
the more digestible 



^/3- 



X 



3/2 



(29) 



TR _ 

For the object masses and semi-major axes used here, this corre- 
sponds to a ratio of 0.758. Therefore, the forcing timescale of ret- 
rograde moons is faster than prograde exomoons, and hence we 
might expect to see more variations in temperature for these orbits. 
Comparing the period of oscillations for cases 1 and 2, we find the 
ratio is indeed approximately 0.758. 

We should also note the amplitude of the oscillations, with the 
prograde displaying an amplitude of around 0.14 K, compared to 
the retrograde amplitude of O.IK. This is again sensible consider- 
ing the orbital dynamics: rp > rp, the prograde moon will spend 
more time at larger and smaller distances from the star, allowing 
the mean temperature to increase and decrease appropriately. The 
retrograde moon will encounter apastron and periastron more fre- 
quently, but spend less time there. This will not allow the mean 
temperature to rise as much, but will facilitate some latitudes to 
increase and decrease their temperature more readily (which gives 
an explanation for the high oscillations seen in the maximum and 
minimum temperatures in Figure|2). 

Case 3 has a significantly lower orbital period (approximately 
3.6 days), compared to the 41 day orbit of the other moons. As the 
diumal and orbital periods are extremely close, and the periastron 
and apastron distances are extremely similar, the effect of epicyclic 
motion is extremely weak, and the climate becomes difficult to dis- 
tinguish from that of an Earth-like planet orbiting with parameters 
(ap,ep). 

The counter-intuitive behaviour of Case 4 becomes clear when 
the alignment of the two orbital angular momentum vectors, Lp 
and Lm are considered. Lp is aligned with the z-axis, and Lm 
is aligned with the x-axis in the negative direction (see Figure |5}. 
At i = 0, the vector product Lp x Lm = L' is aligned with 
the negative y direction. As there are no external torques acting 
on the system, L' maintains this alignment throughout the orbit, 



3.2 Parameter study 

The previous section has shown us that while the global properties 
of the exomoon are relatively insensitive to the orbital architecture, 
the detailed behaviour of the exomoon's climate depends on the 
orbital direction and period. We shall now look at how the global 
properties depend on the parameters more easily determinable by 
exoplanet and exomoon searches: the semi major axes of both ob- 
jects relative to their host, and their eccentricities. To do this, we 
carry out over 3500 separate climate simulations, and classify each 
one as hot, cold, habitable or transient according to the classifica- 
tion system outlined previously. 



3.2.1 A Control - The Circumstellar Planetary Habitable Zone 

To provide grounds for comparison, we ran a separate series of 440 
simulations with an Earth-like planet orbiting a Sun-like star to map 
out the traditional circumstellar HZ according to our classification 
system. The results of this are displayed in Figure |6] Green points 
represent simulations where the planet is classified as habitable; 
purple points where the planet is classified as a transient; and the 
red and blue points indicate planets that are too hot and too cold 
respectively. 

The model does not contain a carbonate-silicate cycle, un- 
like e.g. Iwilliams & Kastind 119971 . which modified the atmo- 
spheric CO2 pressure in accordance with temperature dependent 
weathering rates. Lower temperatures produce lower weathering 
rates, and as a result the atmospheric CO2 (produced by vol- 
canic outgassing) cannot be sequestered. Therefore, cooler plan- 
ets can be expected to have higher atmospheric concentrations of 
CO2, boosting the greenhouse effect and moving the outer edge 
of the HZ to higher semi-major axes than we see in Figure |6] (cf 
iKasting. Whitmlre & Revnolddl 19931) . 

The extension of the HZ at low Cp to semi-major axes as low 
as 0.7 AU is a reflection of our (fairly lenient) criterion for habit- 
ability - namely, that 10% of the planet's surface remains habitable 
over a ten year period, with a standard deviation less than 10% of 
the mean habitable area. As Ep is increased, aj increases quickly, 
producing planets which are habitable on a seasonal basis only. If 
we are to compare to traditional habitability studies, then we should 
infer that the inner edge of the HZ is at approximately 0.8 AU for 
Bp = 0. Equally, many of the transient classifications in this study 
would normally have been considered to be uninhabitable (as many 
of these simulations undergo seasonal periods when the habitable 
surface fraction is close to zero, but this is not sufficient to label 
them as hot or cold planets). We should bear this in mind as we 
consider the habitable zone for exomoons in the following sections. 
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Figure 6. The habitable zone for an Earth-Hke planet around a Sun-like star, as calculated from a LEBM using the classification system outlined in this paper, 
for comparison with the exomoon habitable zones displayed in later figures. Each point represents a simulation run with these parameters, and the colour of 
the point indicates its outcome. Red points produce hot planets with no habitable suiface; blue points produce cold planets with no habitable surface; green 
points represent warm planets with at least ten percent of the surface habitable and low seasonal fluctuations; purple points represent warm moons with high 
seasonal fluctuations. 
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Figure 7. The exomoon habitable zone, as a function of host planet semi-major axis ap and eccentricity ep, for moons with zero orbital eccentricity. Each 
point represents a simulation ran with these parameters, and the colour of the point indicates its outcome. Red points produce hot moons with no habitable 
surface; blue points produce cold moons with no habitable surface; green points represent warm moons with at least ten percent of the surface habitable, and 
low seasonal fluctuations; purple points represent warm moons with high seasonal fluctuations. In each plot, the moon's orbital semi-major axis relative to the 
host planet, a^ is fixed at a value between 0.1 and 0.3 Hill Radii. 
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3.3 Zero Moon Orbital Eccentricity 

We first consider the habitable zone manifold in the case where 
fim = 0, (and consequently the tidal heating rate is zero) and we 
consider several different values of am. This allows us to construct 
Up — Ep slices of the HZ manifold to compare with the typical planet 
HZ shown in the previous section. Figure [7] shows four a^ — Ep 
slices, for four different values of am ■ In each simulation, we fix am 
relative to the Hill Radius at the given value of Op, so the absolute 
value of am will increase with increasing Op. 

Comparing the top left panel of Figure |7] with Figure|6] we can 
immediately see that the habitable zone exists in general at lower 
semi-major axis. This is especially true at high e^: the inner edge 
of the HZ at e = 0.5 is now at 0.84 AU, as compared to 0.89 AU in 
the planetary case. As am is increased, the high Cp component of 
the inner HZ boundary shifts to higher and higher ap. We therefore 
surmise that the effect of frequent eclipses allows the moon to re- 
main cooler, and damps the fluctuations in temperature experienced 
as the planet eclipses the star. 

In all cases, the low Op, low ep habitable tail observed in 
the planetary case disappears in the exomoon case. While frequent 
eclipses can cool the climate and damp oscillations, the habitable 
surface fraction of these moons is generally low (typically close 
to the minimum O.I). Therefore, a relatively small increase in the 
value of G^ due to eclipses is large enough to ensure ctj > 0.1^, 
and as a result be classified as transient moons rather than habitable 
moons. 

Comparing the outer edge of the HZ in the planet and moon 
cases, we can see it takes on a significantly different shape. While 
the planet outer HZ moves steadily outward as Ep is increased, the 
moon outer HZ changes direction as Cp is increased. Initially hab- 
itable at 0.99 AU for Ep — 0, the outer boundary moves inward as 
Bp increases, finding its minimum value of ~ 0.95 AU at Cp = 0.2, 
and then moves outward again. This appears to be due to a combi- 
nation of planet and moon motion relative to the star. The moon's 
epicyclic motion at ep = changes its distance from the star by 
only a small amount, and hence the climate remains stable enough 
to be classified as habitable. As ep is increased, the epicyclic mo- 
tion of the moon is modulated by the motion of the planet as it 
moves between its own periastron and apastron. This extra mod- 
ulation is sufficient to trigger a positive feedback in albedo, and 
a subsequent snowball effect. As Cp increases to large values, the 
strong heating the moon experiences as the planet passes through 
its periastron is sufficient to keep it warm throughout the rest of 
the orbit, and the outer HZ boundary remains at Up ~ 0.99 AU. 
As the eclipse duration increases with increasing am , the outer HZ 
boundary is pushed inwards for larger ep as am is increased. 



3.4 High Moon Orbital Eccentricity 

We now repeat the above analysis, but increase the moon's ec- 
centricity to Bp — 0.1. This is quite a large value to assume for 
moon eccentricity - without eccentricity pumping, it is unlikely that 
moons will possess eccentricities as large as this after tidal evolu- 
tion with their host planet - we present these results as an extreme 
case. Intermediate eccentricities will produce results somewhere 
between the results of this section and the results of the previous 
section. 

Indeed, Figure [8]illustrates how extreme this case is by the ef- 
fect on the climate of the simulations. As tidal heating increases as 
a,„ is decreased, it is clear that tidal heating dominates the proper- 



ties of highly eccentric moons for am < 0.2 Hill Radii (top panels), 
and all simulations are classified as hot. 

As am is increased, a habitable zone appears. This zone is 
much narrower than seen previously, and is restricted to large val- 
ues of ap, as tidal heating allows for warmer solutions when inso- 
lation is weak. In the case where am = 0.21 Hill Radii (bottom left 
panel of Figure[8), the HZ exists at ap > 0.95 AU, and potentially 
extends beyond I.I AU at high ep. The HZ is typically around 0.1 
AU in width, although at low ep it is as narrow as 0.03 AU. 

When am is increased to 0.26 Hill Radii, the HZ begins to 
more closely resemble the HZ seen for e™ = moons. While 
still around 0. 1 AU in width, it now displays the knee in the outer 
HZ boundary seen previously, and it is centred on similar values 
of ttp as before. At this value of am, the tidal heating has become 
small compared to the traditional heating from insolation. This be- 
ing the case, tidal heating still makes its presence felt at low ap. 



where simulations that would have been habitable for Cr, 



are 



now seasonally habitable due to the extra forcing from tidal heating 
producing stronger climate oscillations. 



4 DISCUSSION 

4.1 Limitations of the Model 

As with all studies using ID LEBMs, the nature of the assumptions 
made limit the applicability of the results obtained in this analy- 
sis. The models are insensitive to the longer term climate varia- 
tions produced by oceanic circulation, as well as the compensatory 
effects of the carbonate-silicate cycle which can halt the albedo- 
driven snowball effect. These effects, with timescales ranging be- 
tween a few thousand to a million years JWilliams & Kastinal997h 
are much longer than the orbital periods of both the moon and the 
host planet, and as such are of less consequence (in the short term) 
than the dynamical forcing due to eclipses and tidal heating. 

Assuming a value for the diffusion constant, D limits the 
application of LEBMs to exomoons in a different sense. D is 
calibrated such that a fiducial Earth-Sun climate model produces 
the corre ct latitudinal temperature distrib ution as obtained from 
real data dSpiegel. Menou & Scharij|2008l) . Therefore, D contains 
a great deal of hidden information about the atmospheric and ther- 
modynamic properties of the Earth. Modifying these properties, 
e.g. the atmospheric press ure, has significant consequences for the 
habita bility of the object JWilliams & Kastinel[l997l : IVladilo et alj 
l20I3h and it is less clear how D ought to be modified for frozen 
moons such as Europa, or moons with fundamentally different 
chemistry such as Titan, without more data on their own latitudinal 
temperature distributions. Even with this data, a full understanding 
of how tidal heating affects the temperature balance at any latitude 
is essential for the correct calibration of D. Equally, a more accu- 
rate implementation of tidal heating in the LEBM is vital for future 
studies of exomoons. 

Tidal forces are expected to affect land and ocean according 
to their elastic rigidities. While we have assumed a fixed ocean 
fraction of 0.7 throughout this analysis, the LEBM does not con- 
tain information as to how these oceans are distributed across the 
surface. Instead, the ocean and land components are assumed to 
be uniformly mixed over all latitudes. Future calculations of tidal 
heating will require the LEBM to be more specific about the dis- 
tribution of landmasses. Ocean fraction may therefore become an 
important parameter in future studies of the exomoon HZ mani- 
fold (which should also investigate the effects of changing other 
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Figure 8. The exomoon habitable zone, as a function of host planet semi-major axis ap and eccentricity ep, for moons with eccentricity em = 0.1. Each 
point represents a simulation ran with these parameters, and the colour of the point indicates its outcome. Red points produce hot moons with no habitable 
surface; blue points produce cold moons with no habitable surface; green points represent warm moons with at least ten percent of the suiface habitable, and 
low seasonal fluctuations; purple points represent warm moons with high seasonal fluctuations. In each plot, the moon's orbital semi-major axis relative to the 
host planet, a^ is fixed at a value between 0.1 and 0.3 Hill Radii. 



parameters such as obliquity). Besides its effect on tidal heating, the 
ocean fraction also affects the thermal relaxation timescale of the 
moon - reducing the ocean fractio n leaves the moon more suscep- 
tible to larger climate oscillati ons ( ISpiegel. Menou & Scharn2008l : 
lAbeetalJl201ll : lForganll2012h . 

Tied to this problem of limited surface data is the limited di- 
mension of the model itself - ignoring longitudinal information re- 
quires us to assume the moon is rapidly rotating relative to the inso- 
lation source. This will remain true for exomoons which are tidally 
locked to the planet - as they will still rotate with respect to the pri- 
mary insolation source, the star - but the more complicated nature 
of the lunar day compared to a planetary day may mean that the 
rapid rotation assumption is satisfied more weakly at some points 
in the orbit compared to others. 

We have also ignored one important heat source - illumina- 
tion of the moon by the planet. iHeller & BamesI ( 120130 show cal- 
culations for illumination from a tidally locked planet, incorporat- 
ing the planet's thermal radiation and reflected starlight, with re- 
flected starlight typically dominating unless the planet's albedo is 
very low. Adding the reflected starlight alone can increase the outer 
HZ boundary by a few percent - if this was incorporated into the 
LEBM, it could act to prevent the "knee" in the outer boundary at 
moderate values of Cr,. 



Finally, we assume orbital parameters for all objects involved, 
and we do not evolve these parameters under the presence of any 
gravitational field. Incorporating the LEBM component into a more 
involved com putation of the tida l evolu t ion of the star-planet-moon 
system (e.g. iLaskar & Robutell 1 19931 ; ISasaki. Barnes & O'BrienI 
l2012b would allow an investigation of a unique set of Milankovitch- 
esque cycles produced by the moon's orbital evolution (see e.g. 
Spiegel et al. 2010), with an accurate computation of the moon's 
obliquity evolution and tidal heating as a bonus. 



4.2 Implications for Observations 

Our simulations suggest that determining the orbital eccentricity 
of any exomoons detected in or near the conventional HZ will be 
crucial in assessing whether the moon itself is also habitable, much 
more so in fact than determining the host planet's eccentricity. Exo- 
moons with even a moderate amount of eccentricity of ~0. 1 lead to 
a more limited habitable-zone with respect to the other system pa- 
rameters as a result of tidal heating dominating. One may consider 
that the HZ extends outwards in ap for such scenarios but longer 
period planets are considerably rarer in the known transiting exo- 
planet catalogue, due to both geometric bias and lower composite 
signal-to-noise. 
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The orbital eccentricity of a moon will affect the position and 
duration of moon transits in the light curve (Kipping 2011) and 
is, in principle, an observable quantity. It is therefore quite reason- 
able that constraints on exomoon eccentricity could be determined 
in the case of a confirmed detection. It is also worth noting that 
exomoons are expected to rapidly tidally circularize even if start- 
ing with a large i nitial e ccentricity due t o a capture, for example 
IPorter & Grundvl ilOl lb : Iwilliamsl ( l2013l) . 

To maintain a large exomoon eccentricity over Gyr would 
likely require forcing due to resonant massive satellites, but the 
occurence of multiple captures of Earth-mass moons around a sin- 
gle planet is surely less than that of single captures of Earth-mass 
moons. 

We also find that the sense of orbital motion affects the cli- 
mate of m oons and in prin ciple this may also be determined obser- 
vationally (Kipping 2009b). Ho wever, in general such a determina- 
tion is highly challenging with lKipping et alj (201J) finding even 
optimstic synthetic data is marginal for uniquely determining this 
to high-confidence using Kepler. Of course, once a confirmed sig- 
nal has been found, follow-up with larger instrumentation such as 
JWST would greatly increase our ability to make this assesment. 

Finally, in this work we have established an LEBM capable 
of modelling hundreds of realizations of exomoon configurations. 
Due to the likely low signal-to-noise of any future exomoon de- 
tection one should expect fuzzy and co mplex parameter po steri- 
ors, as demonstrated in synthetic tests bv lKipping et aljj2012h . We 
therefore suggest coupling an LEBM, such as that presented here, 
directly with representative realizations drawn from the j oint pos- 
terior distribution of the relevant parameters. This would enable us 
to compute a marginalised distribution for the habitability of an ex- 
omoon (^). 



5 CONCLUSIONS 

In this work, we have used a one dimensional latitudinal energy 
balance model (LEBM) to investigate the evolution of climate on 
an Earth-like exomoon. We have incorporated the effects of stellar 
insolation, tidal heating, eclipses of the star by the planet, atmo- 
spheric cooling and heat diffusion, and investigated how the dy- 
namical circumstances of exomoons affect the resulting climate. 

We simulated four test cases to study the exomoon climates 
in detail, finding that exomoons in orbits retrograde to their host 
planet's orbit experience stronger climate oscillations than their 
prograde equivalents, and are around O.IK warmer. In the case of 
moons with sufficiently short orbital periods, we find that the fi- 
nite thermal relaxation time of the atmosphere allows it to act as 
a buffer against the otherwise strong rapid temperature oscillations 
produced by frequent eclipses. 

We then went on to carry out an extensive parameter study of 
the four dimensional space constituted by the planet's semi-major 
axis and eccentrictity (Op, e^) and the moon's semi-major axis and 
eccentricity relative to the planet (a™ , e^). Comparing to the habit- 
able zone of an Earth-like planet orbiting the same star, we find that 
if the exomoon's orbit is circular, the exomoon habitable manifold 
tends to exist at lower ap thanks to the cooling effect of eclipses, 
and the inner edge of the habitable zone can exist at lower Op for 
higher ep, again as eclipses keep the moon sufficiently cool. If the 
exomoon's orbit is highly eccentric, the heating budget is domi- 
nated by tidal effects, and the habitable manifold is much narrower 
in flp, and exists at greater distance from the star. 

The simplicity of LEBMs, plus their (albeit limited) ability to 



be augmented by new heat sources and sinks, makes them a useful 
tool for studying Earth-like exomoon climates. As they are com- 
putationally inexpensive, it is straightforward to run them multiple 
times, whether it is to map out their behaviour in a set parameter 
space as done here, or to assess the habitability of detected exo- 
moons, by using realisations of the joint posterior distribution de- 
rived from observations. We are keen to hone and utilise this mod- 
elling technique for both of these uses. 
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